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Abstract: In this paper, we compare the spatially adaptive 
iterative filtering (SAIF) approach with Weighted Encoding 
with Sparse Nonlocal Regularization (WESNR) to maintain 
the denoising strength locally for any spatial domain method. 
These approaches has ability of filtering local image content 
iteratively using the given base filter, and the type of iteration 
and the iteration number are automatically optimized with 
respect to estimated risk (i.e., mean squared error). 
Experimental result shows that approx. 7% enhancement in 
the SNR of SAIF method as compared to the WESNR method. 

Keywords: Image denoising, pixel aggregation, risk 
estimator, spatial domain filter. 

I. Introduction 

SINCE noise is a random phenomenon which occur basically all 
modern imaging systems in the process of transmission and 
reception, and denoising is used to restoration of fundamental 
image. There has been invented a lot number of denoising 
algorithms, and in general they are categories in two main parts: 
first one transform domain methods and second one spatial 
domain methods. Transform domain methods are used as a 
image is represented as a combination of certain transform 
basis function, by which signal to noise -ratio can be calculate 
and used the appropriate shrink for respective transform 
coefficient. 

In Non-local mean (NLM) as we know it is data 
dependent filter, each pixel is calculated as the weighted 
average of all its similar pixels in the image, and by using 
similarity between them determined their weights. Zhang et al. 
[15] clustered the similar patches into a matrix and applied 
principal component analysis (PCA) to eliminate AWGN. The 
so-called LPG-PCA algorithm provides very good edge 
preservation performance. But in current years, In image 
restoration the sparse representation and dictionary learning 
based methods have been attracting significant attention. The 
combine use of sparse representation and nonlocal self- 
similarity regularization has led to idea of AWGN removal [26] . 
Another successful method called non-local means (NLM) 
extends the bilateral filter by replacing point-wise photometric 
distance with patch distances, which is more robust to noise [6]. 
In general, for spatial domain method the main problem is to 
determine the denoising strength. Because in this method 
denoising strength mainly depends on tuning parameter .A low 
value of controlling parameter would maintain high-frequency 
signal but will do minute denoising (estimation variance) in the 



image. And larger value of controlling parameter would quash 
more noise and also remove some important image information, 
and finally yield over-smoothed image. Another option for this 
controlling parameter selection is iterative filtering. By using 
this technique, either a filter estimated by wrong tuning 
parameter, but by using this filter method several numbers of 
times, can still get a well estimated output then iteration 
number then becomes another tuning parameter that needs to be 
carefully treated .In the spatial domain method, denoising 
strength is automatically adjusted according to local SNR. In 
[7] Milan far explains that a spatial domain denoising process 
can be equivalent as a transform domain filter, where the 
orthogonal basis elements are the eigenvectors of a symmetric 
and positive definite matrix determined by the filter; and the 
shrinkage coefficients are the corresponding eigenvalues 
ranging in [0, 1]. For NLM filters the eigenvectors 
corresponding to the dominant eigenvalues could well represent 
latent image contents. According to this idea, we suggest a 
spatially adapted iterative filtering (SAIF) strategy capable of 
controlling the denoising strength locally for any given spatial 
domain method. The proposed method iteratively filters local 
image patches, and the iteration method and iteration number 
are automatically optimized with respect to local MSE, which 
is estimated from the given image. To find out the MSE for 
each patch, we discuss plug-in risk estimator. While [9] also 
uses SURE to optimize the NLM kernel parameters, we 
illustrate that (1) the plug-in estimator can be superior for the 
same task, and (2) the adaptation approach can be extended to 
be spatially varying. This review paper is arranged as in Section 
II will provide some background, especially [7]'s proper 
analysis of spatial domain filters. Section III reviews two 
iterative methods to regulate the controlling (smoothing) 
strength for the filters in section III . Section IV describes SAIF 
strategy in detail. Experimental results are given in Section V to 
show the performance of the SAIF strategy using several 
leading filters. Finally we conclude this paper in Section VI. 

II. Background 

Assume the typical measurement model for the denoising 
difficulties: ?i = z i + e i fffr i = 1 - (1) 

where = z(x^ ) is the positioned image at 

position X{ = [^,1,:^, 2] r , y- t is the noisy pixel value, and Sj 
denotes zeromean white noise with variance u*. The problem 
of denoising is to recover the set of underlying 
samples z = \z v . . . , z n ] T . The representation of vector model 
in the measurement model is: 

y = z + 9 (2) 
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As explained in [5], [7] most spatial domain filters can be 
represented through the following non-parametric restoration 
framework: 

z { = argminl^L^ - y^] 2 K(x if xj f y if yj') (3) 

where stands the estimated pixel at position and the 
weight (or kernel) function K(-) measures the similarity 
between the samples y\ positioned at X; and y.- positioned at Xy 

In general the widely used kernel function is the Bilateral (BL) 
filter [4], that smooth images by the help of a nonlinear 
combination of nearby image values. The method combines 
pixel values based on both their geometric closeness and their 
photometric similarity. This kernel can be expressed in a 
separable fashion as follows: 



tfij = exp 1 ~, + 



(4) 



in which h x and hy are controlling (smoothing) parameters. The 

closely look like the bilateral filter except that the photometric 
similarity is collected in a patch wise manner: 

-hi-yj 



Kij = gxp 



r- : : 



+ 



r * 



(5) 



Where yt patch positioned at y ? i and yj patch positioned at yj . 
Only in theory (not in practical way,) the NLM kernel has the 
patch- wise photometric distance (h^ — >oo). 
In common way, preceding restoration algorithms are based on 
the same framework (3) in which some data-adaptive kernels 
are assigned to each pixel contributing to the filtering. 
Minimizing equation (3) gives a normalized weighted averaging 
process: 

Zj = wjy (7) 

where the weight vector Wj is 



W £ = 



(8) 



L. 



Collecting the weight vectors together, the filtering process is 
proceed for all the sample pixels simultaneously by a 
multiplication form of matrix-vector 



z = 



T-i 



T 



LU"£ J 



y = Wy 



(9) 



Where z denote estimate signal and W denote the filter matrix. 
W is a positive row- stochastic matrix (sum of every row in 
matrix will be 1). This matrix will not be symmetric matrix but 
consist of positive real eigenvalue [7]. Even W is not a 
symmetric matrix in general, it can be closely estimated with a 
symmetric positive definite matrix [12]. The symmetrized W 
must also stay row- stochastic, that means we get a symmetric 
positive definite matrix which is doubly (i.e., row- and column-) 
stochastic. Thus we can estimate eigen -decomposition of 
symmetric W as follows: 

W = VSV T (10) 
In this equation S = diag[A LJ ... ... , A n ] consist of the eigenvalues 

in descending order 0 < A L <, . . . ,< A n = 1, and V is an 

orthogonal matrix V = [v v f v n ] consist of the corresponding 

eigenvectors of W in its columns. Because V is orthogonal, its 
columns specify a set of basis functions, therefore the filtering 
process can be expressed as: 
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$ = Wy = VSV T y (11) 
Here by using the eigenvectors of W the input data y is 
convert into the domain spanned; then, each coefficient is 
modify (scaled) by the factor A - t ; and an inverse transformation 
is done in final step, to get the output. By doing the above 
analysis we find that the denoising strength for every basis of a 
given filter is thus controlled by the scaled factor {Aj}. In the 
next sections we will elaborate this concept more broadly. 

III. Iterative Filtering Methods 

Main shrinkage procedures based on various spatial domain 
filters have been discussed in [7], where statistical analysis 
shows that the optimal filter with respect to MSE is the local 

Wiener filter wither = it, where snr; indicate 



ith channel signal-to-noise ratio. In, the local Wiener filter 
require proper knowledge of the local signal to noise (SNR) of 
every basis channel, that does not directly accessible in general. 
In denoising schemes such as [1] and [3] the Wiener shrinkage 
criterion works centred on a pilot estimate of the latent image. 
Still, the Wiener filter's performance firmly believes on 
accuracy of this estimate. Iterative filtering can be a trustworthy 
substitute for reducing sensitivity of the basis shrinkage to the 
calculated local SNR. Then, the iteration number will be the 
only parameter to be locally optimized. To approach the locally 
optimal filter performance in a stable way, we propose the use 
of two iterative local operators; namely diffusion and boosting. 
In [13] we have shown that performance of any type of kernel 
could be enhanced by iterative diffusion which gradually 
removes the noise in each iteration. Yet, diffusion filtering also 
takes away latent details from the underlying signal. On the 
other hand, iterative boosting is a mechanism to preserve these 
lost details of the signal. By using the two iterative filtering 
methods, we can avoid either over-smoothing or under- 
smoothing due to incorrect parameter settings. In another way, 
these two methods deliver a way to start with any filter, and 
maintain the values of shrinkage factors {A{} to find a better 
and firm approximation of the Wiener filter. Further we have 
explained the two methods, separately. 
A. Diffusion 

The notion of diffusion in image filtering was originally 
followed by the basic idea of heat propagation and explained 
using a partial differential equation, this method adopt the 
discrete version of it, that is easily represented by repeated 
applications of the same filter as described in [7] 

Each application of W can be interpreted as one step of 
anisotropic diffusion with the filter W.By taking large value of 
iteration number of k will eliminate the noise and high 
frequency component of noise, over smooth the image at the 
same time and small number of iteration k maintain the 
underlying structure. Minimization of MSE (or more accurately 
an estimate of it) stands the proper time to stop filtering, that 
will avoid becoming under- or over- smoothing of image. 
Here W is symmetric, the filter in the iterative model (12) can 
be decomposed as: 

w k _ vs k v T (13) 
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Here = diagP.*, . . . It is important that in spite of the 
general prediction of k as a discrete step, by the realization 
spectral decomposition of makes it possible to substitute k 
with any positive real number. The hidden image z would be in 
the column space of V as b = V T z, where b = [bl, b2, . . . , bn] 1 
, and {bt} stands the projected signal energy over all the 
eigenvectors. As shown in [7] the iterative estimator z- A = W^y 
has the following squared bias: 

\\b iaSk ||* = ||(J - H^)zP = EF =1 (1 - X^bf (14) 

Correspondingly, the estimator's variance is: 

var(2 t )=tr(cov(^))=cr : X^f (15) 
Finally, the MSE is given by 

MSE k =|| bias k || 2 +var (z ) = S? = L ( 1 - A*) 2 bf + a 2 Xf ( 1 6) 
As the iteration number k grows, the bias term increases, but the 
variance decays to the constant value of u \ Of course, bthis 
expression for the MSE is not practically useful yet, since the 

n. 

coefficients {b\ } are not known. Later we describe a way to 
estimate the MSE in practice. But first, let us introduce the 
second iterative mechanism which we will employ. Boosting is 
discussed in the following and as we will see, its behaviour is 
quite different from the diffusion filtering. 
B. Boosting 

Meanwhile diffusion filtering is used in more practice but in 
low SNR case diffusion filtering method is fail. And in that 
case, diffusion always eliminates some part of the noise and 
signal, concurrently. This problem is short out by boosting that 
recycles the removed components of signal from the residuals, 
in each iteration. Defining the residuals as the difference 
between the estimated signal and the noisy signal: 
r k = y — the iterated estimate can be expressed as 

= (i- (i- wo ft+1 )y ( 17 ) 

where z$= Wy. And as k increases, the estimate backs to the 
noisy signal y. basically the boosting filter has different 
behavior than the diffusion filter where after each iteration the 
calculated signal gets close to a constant. After k iteration the 
squared magnitude of the bias is 

\\bias k \\ z = ll(/- HO ft+1 z|| a =Sf =1 (l -X[) 2k + 2 bf (18) 
Estimators' variance is, 

Var(z k )=tr(cov(f k ))=a 2 £? =1 [1 - (1 - A £ ) ft+1 ] 2 (19) 
the overall MSE is 

MSE k = S? =1 (l - li) ait+3 &f + a\l - (1 - JL £ }* +1 ) a (20) 
As k increases, the variance will increase and bias term will 
decrease. Highlighting the nature of the diffusion iteration, 
analyse that boosting replace diffusion when diffusion unable to 
improve the filtering performance. The support of this task is 
that we can automatically optimize the type and number of 
iterations locally to boost the performance of a given base 
filtering method at the same time. 
Algorithm of WESNR 

Here, the dictionary <p is given, assume. It is an important issue 
selection of dictionary to the sparse coding and reconstruction 
of a signal. In specific, from natural image, learning dictionaries 
patches has given favourable results in image restoration [17, 
18]. In this algorithm, the dictionary <p is pre-learned from clean 
natural images, and the pixels corrupted by Impulsive Noise 
(IN) will have big coding residuals. Consequently, the coding 
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residual can be used to guide the setting of weight Wa , and 
should have inverse relation with the strength of &i . In order 
to make the weighted encoding stable and easy to control, we 
select Wii should lie [0; 1]. So, relation between and &i is 

W£ = expi-a*fi; (21)' 
Here should be decreasing rate with respect to &i.To 
maintain this relation a will be positive constant By Eq. (21), 
the pixels corrupted by IN will be adaptively assigned with 
lower weights to reduce their impact in the process of encoding. 
Let V be a diagonal matrix and initialize it first as an identity 
matrix, and then in the (k + l)th iteration, each element of V is 
restructured as 

1$ *i = !/[(«* - ft ) 3 + ^] f £ (22) 

where e is a scalar and a- ^ is the E th element of coding vector 
a in the k th iteration. Then we update cl as 
gk + i = ^T W(f +yto+tiy\ tp Twy - 9 T H/^ + /j (23) 

By iteratively updating V and a, the desired a can be efficiently 
obtained. A number of patches (size: 7X7) are extracted from 
the five images and they are clustered into 200 clusters by using 
the K-means clustering algorithm. For each cluster, a compact 
local PCA dictionary is learned. Meanwhile, the centroid of 
each cluster is calculated. For a given image patch, the 
Euclidian distance between it and the centroid of each cluster is 
computed, and the PCA dictionary associated with its closest 
cluster is chosen to encode the given patch. Note that since the 
selected dictionary, denoted by^i, is orthogonal, the^r for 
patch X{ can be simply computed as u = <pj 
First the dictionary <p is adaptively calculated for a given patch, 
the WESNR model can be solved by iteratively updating W 
and a. The updating of W depends on the coding residual e. In 
the literature of mixed AWGN and SPIN noise removal [19]- 
[22], AMF [2] is widely used to detect SPIN. In order to make a 
fair comparison with them, in the case of AWGN+SPIN noise 
removal, we apply AMF to y to obtain an initialized image 
-f(0], and then initialize e as: 

g m = y - x W (24) 

In the case of AWGN+RVIN+SPIN noise removal, AMF 
cannot be applied to y to initialize x. We initialize e as 

■ TO = y-f*j»i (25) 

where p y is the mean value of all pixels in y and 1 is a column 

vector whose elements are all 1. In other words, we simply use 
the mean value of y to initialize x. Then the initial coding 
residual can be roughly computed. This simple initialization 
strategy works very well in all our experiments. With the 
initialized coding residual W can be initialized by Eq. 
(21). The main procedures of the WESNR based mixed noise 
removal algorithm are brief in Algorithm 1. In our algorithm, 
we set t = \\^a^ - * ff M||/||* ff m|| < T as the 

termination condition. 

The code for the WESNR algorithm as: 

Algorithm 1: WESNR method to eliminate mixed noise 

Input: Dictionary^, noisy image y ; 

Initialize e by Eq. (24) (or Eq. (25)) again 

Initialization of u to 0 

i 

Initialization of W by Eq. (21); 
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Output: Denoisey image x. 
Loop: iterate on k = 1,2, . . .,K; 

1. Find out by Eq. (23); 

2. Find out cc^ = tptx^ and accordingly change the nonlocal 
coding vector u ; 

3. Find out the residuals ^ = y — x ^ ; 

4. Compute the weights value W by using Eq. (21); 
End 

Output will be denoisey image x = tpa^. 

IV. Modified Saif Method 

Based on the analysis from Section III we propose an image 
denoising strategy which, given any filter using the framework 
(3), can boost its performance by utilizing its spatially adapted 
transform and by employing an optimized iteration method. 
This iterative filtering is implemented patch wise, so that it is 
capable of automatically adjusting the local smoothing strength 
according to local SNR. 
A. Optimal Iteration Estimation 

Let us assume a patch y and corresponding filter matrix W, here 
main aim is to select appropriate iteration method and number 
of iteration that will minimize the MSE. In brief, the optimal 
stopping time k for each iteration method is given as: 

k 

The method to calculate an unbiased estimate of MSE is the 
SURE [8]. The best alternative, here we propound, is the plug-in 
risk estimator, that is biased and works based on an estimate of 
the local SNR. First, we estimate, eigenvalues and eigenvectors 
of the filter from a pre-filtered patch z, obtained with the help 
of the base filter by setting some arbitrary parameter. We have 
expressed as: 

W(z') = VS V J (22) 
In spite of the earlier prediction of k as a discrete step, in 
practical way the spectral decomposition of W indicate that, k 
can be any positive real number. To be more 

definite, W k = VS k V T , with S k = diag[ a\ 1*] where k is 

positive real number. In practical implementation, the filter can 
be suggested with modified eigenvalues for any k > 0. This will 
improve the performance .In consequence; a real- valued k 
automatically takes and smoothly adjusts the local bandwidth of 
the filter. While decreasing the iteration number k can be 
predicted as smaller tuning parameter fr v for NLM kernel, larger 

the value k is equivalent to a wider kernel. Further, we explain 
the two risk estimators and find out that the plug-in estimator 
will give better performance in comparison to the SURE 
estimator for estimated local SNR. 

1 ) Plug-In Risk Estimator: It is explain in Algorithm 2. In this 
technique, risk estimators for boosting and diffusion are 
calculated based on the pre filtered patch z, computed using the 
base filter with arbitrary parameters. More clearly, for the 
estimation of the signal coefficients as follow: 

B=V T z (23) 
This equation's contribution can be predicted as equipping the 
risk estimator with some prior idea of the local SNR of the 
image. The estimated signal coefficients allow us to use (24) 
and (25) to find out MSE-^ in each patch: Diffusion Plug-in Risk 
Estimator is given as: 



Plug - in( = £[^(1 - + a*J% k (24) 

Boosting Plug-in Risk Estimator is given as: 

ping - m-° = u =1 (i - 4> 2ft+2 s? + ff 3 d - (i - ^) i+i y 

(25) 

In each patch, minimum values of Plug — and Plug — in^ 
as a function of k are calculated and then compared, and the 
iteration type with the least risk is selected. Since the optimal 
iteration number fccan be any real positive value, in the 

application of the diffusion filter, W is substitute by VS* 1 V 1 y 

s s ^ 
where S* 1 = diag A£ ]. This has been similarly 

shown for the boosting filter in Algorithm 1 . Next, for the sake 

of comparison, the SURE estimator is discussed. 

Algorithm 2: Plug-in Risk Estimator 

Input: Noisy Patch: y, Pre-filtered Patch: z, Patch 

Filter: W 

Output: Denoised Patch: z 



1. 

2. 
3. 

4. 
5. 

6. 
7. 

8. 



Eigen-decomposition of the filter W(z ) = l^5^ T 
B = V* z t= Compute the signal coefficients; 
Plug - irtft , Plug - in® <= Compute the estimate 
risks; 

if min{F^ - in£} < [Plug - ing} 
k.= arg m * n Plug - wt^<= Diffusion optimal 
iteration number; 

z = VS^V 7 y <= Diffusion patch denoising; 
else 

■■ _ argmin 



Plug — mf <= Boosting optimal 



iteration number; 

z = V{l - (1 - S} £ +1 )t^y t= Boosting patch 

denoising; 
End 

V. Experimental Results 

In this section we evaluate performance of method for 
denoising various images .we first compare results with ones 
from the standard kernels: NLM and Bilateral. 



9. 



10. 
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Fig. 10. Comparison of denoising performance on noisy Parrot image corrupted 
by AWGN of a = 25. (a) Original image, (b) Noisy input, (c) WESNR [6]. (d) 
SAIF (NLM). 

TABLE-I 

The PSNR Comparison between WESNR and SAIF Methods 



Method 


PSNR 


Time (Sec.) 


Noisy Image 


20.15 




WESNR 


27.18 


23.37 


SAIF 


28.87 


90.22 



We also test stability of SAIF when an arbitrary tuning 
parameter is used. In all cases the estimators show a promising 
improvement over the standard kernels. We will show that the 
resulting SAIF-ly improved filters are comparable, in terms of 
MSE (PSNR) and SSIM [14], to advanced denoising methods, 
and in many cases visually better. 



PSNR 




■ PSNR 



Noisy Image WESNR SAIF 



As shown, the diffusion gets better most of the flat and smooth 
patches, boosting cares the texture and more complex ones. It is 
pointed that applying overlapped patches has the advantage of 
computing multiple estimates for each pixel from the both 
iteration types. In another way, in the aggregation process, some 
pixels may be from diffusion in one patch, whereas others may 
be from boosting in another overlapping patch. This is 
intentionally helpful for pixels at the border of smooth and 
texture regions. 

V. Conclusion 

For any spatial domain filter, we can boost its performance 
highly developed by employing optimized iteration methods. 
Patch-wise implementation is followed in iterative filtering. 
Possessed with diffusion and boosting as two complementary 
iteration techniques, the optimum local filter filtered each patch. 
More precisely, by using the best iteration number and method 
that minimizes MSE in each patch, SAIF is capable of 
automatically adjusting the local smoothing strength according 
to local SNR. The weighted encoding (WESNR) technique gets 
rid of the Gaussian noise and impulse noise at a time. But the 
final result in terms of PSNR the SAIF gives 7% better 
performance. 
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